Isopy in action!¶
Testing Isopy's three functions using data from Argo float products
In [1]:
import numpy as np
import xarray as xr
from matplotlib import pyplot as plt
%matplotlib inline
In [2]:
url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/Argo_Products/monthly_mean/Gridded_monthly_mean'
ds = xr.open_dataset(url)
ds
Out[2]:
In [3]:
# %load ../isopy/isopy.py
import numpy as np
import xarray as xr
from scipy import interpolate
def compute_iso(data, iso, lev, d_or_t):
""" compute the depth of some isopycnal/isotherm.
Parameters:
----------------
data: the input variable (Temperature, Density),
array_like (3D), shape (N, M, L).
iso: the objective isosurface layer,
float.
lev: the vertical depth,
array_like (1D), shape (N).
d_or_t: input d or t for density or temperature (lowercase)
string
Returns:
----------------
var_iso: the output variable,
array_like (2D), shape (M, L).
"""
size = data.shape
L = size[2]
M = size[1]
N = size[0] #this should be the same length as lev
#check that the input data is in the for depth lat lon (or depth lon lat)
if len(lev) != N:
print('Error: either data input variable does not have depth as dimension 0, or lev variable is not same length as depth')
else:
if d_or_t in {'d','t'} :
var_iso = np.zeros((M, L)) # define the output var
var_iso[:, :] = np.nan # NaN fill to avoid any later computation errors with zeros
for i in np.arange(L): #loop through dimension 1
for j in np.arange(M): #loop through dimension 2
data_prof = data[:, j, i] #select one profile
if d_or_t == 'd':
id1 = np.where(data_prof < iso)
if np.size(id1) > 0 and np.size(id1) < len(z):
var1 = data_prof[id1[0][-1]]
var2 = data_prof[id1[0][-1] + 1]
if var1 < var2:
func = interpolate.interp1d([var1, var2], [lev[id1[0][-1]], lev[id1[0][-1] + 1]])
var_iso[j, i] = func(iso)
else:
id1 = np.where(data_prof > iso)
if np.size(id1) > 0 and np.size(id1) < len(z):
var1 = data_prof[id1[0][-1]]
var2 = data_prof[id1[0][-1] + 1]
if var1 > var2:
func = interpolate.interp1d([var1, var2], [lev[id1[0][-1]], lev[id1[0][-1] + 1]])
var_iso[j, i] = func(iso)
else:
print('Is this density or temperature?')
return var_iso
def iso_val(density_data, ts_data, iso):
""" compute temperature or salinity along some isopycnal/isotherm.
Parameters:
----------------
density_data: density input,
array_like (3D), shape (N, M, L).
ts_data: temp or salt input,
array_like (3D), shape (N, M, L).
iso: the objective isosurface layer,
float.
Returns:
----------------
ts_iso: the output variable,
array_like (2D), shape (M, L).
"""
size = density_data.shape
L = size[2]
M = size[1]
N = size[0]
ts_iso = np.zeros((M, L)) # define the output var
ts_iso[:, :] = np.nan # NaN fill to avoid any later computation errors with zeros
dens = xr.DataArray(density_data)
var = xr.DataArray(ts_data)
a = dens.where(dens.min(axis=0)<iso)
b = var.where(dens.min(axis=0)<iso)
for j in range(M):
for k in range(L):
x = a[:,j,k]
y = b[:,j,k]
f = interpolate.interp1d(x.values,y.values,'linear')
ts_iso[j,k] = f(iso)
return ts_iso
def iso_average(data, var_iso, lev):
""" compute vertical average above the objective isopycnal layer
Parameters:
----------------
data: the input variable (temperature, salinity, U, V),
array_like (3D), shape (N, M, L).
var_iso: the isopycnal layer,
array_like (2D), shape (M, L).
lev: the vertical depth,
array_like (1D), shape (N).
Returns:
----------------
var_ave: the output variable,
array_like (2D), shape (M, L).
"""
size = data.shape
L = size[2]
M = size[1]
N = size[0]
data_new = (data[:-1, :, :] + data[1:, :, :]) / 2
dz = lev[1:] - lev[:-1]
var_out = np.zeros((M, L))
var_out[:, :] = np.NaN
for i in np.arange(L): #loop through dimension 1
for j in np.arange(M): #loop through dimension 2
var1 = data_new[:, j, i] * dz
idx = find_nearest(lev, dep_iso[j, i])
if idx > 0:
var_out[j, i] = np.nansum(var1[:idx]) / np.sum(dz[:idx])
return var_out
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
Test Cell for Func 1: compute_iso¶
In [4]:
%%time
# read density, temp, salinity
pden = ds.pden #.chunk(chunks={'lat' : 30, 'lon' : 30})
temp = ds.temp
sal = ds.salt
z = pden.lev.values
iso = [27]
for t in np.arange(1):
# isopycnal depth
data0 = pden[90, :, :].values
dep_iso = compute_iso(data0, iso, z, 'd')
# now, try isothermal depth
data1 = temp[90, :, :].values
temp_iso = compute_iso(data1, iso, z, 't')
Now, Let's plot it in Cartopy¶
In [5]:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def make_map(projection=ccrs.PlateCarree(), figsize=(20, 15)):
fig, ax = plt.subplots(
figsize=figsize,
subplot_kw={"projection": projection})
return fig, ax
In [55]:
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
lons, lats = np.meshgrid(lon, lat)
h = ax.contourf(lons, lats, dep_iso, trans=ccrs.PlateCarree())
ax.contour(lons lats, dep_iso, trans=ccrs.PlateCarree(), colors='w')
plt.colorbar(h, orientation='vertical', fraction=0.02).set_label(label='ËšC', size=15)
ax.set_title('Depth field of the \nabove 90 isopycnal mark', size=20);
Now, let's try hvplot¶
In [7]:
# change numpy to Xarray
ds_dep_iso = xr.DataArray(dep_iso, coords={'lat': lat, 'lon': lon}, dims=['lat', 'lon'])
ds_dep_iso
Out[7]:
In [9]:
import hvplot.xarray
import geoviews as gv
import holoviews as hv
from holoviews import opts
In [54]:
proj = crs.PlateCarree()
ds_dep_iso.hvplot.quadmesh(
'lon', 'lat', projection=proj, project=True, global_extent=True,
width=1000, height=540, cmap='magma_r', rasterize=True, dynamic=False
) * gv.feature.coastline * gv.feature.land(color='red')
# ds_dep_iso.hvplot.contour('lon', 'lat', projection=proj, project=True, global_extent=True,
# width=1000, height=540, rasterize=True, dynamic=False, line_width=5
# ) * gv.feature.coastline
Out[54]:
Test Cell for Func 2: iso_val¶
In [12]:
%%time
isotemp = iso_val(pden[90, :, :].values, temp[90, :, :].values, iso)
isosal = iso_val(pden[90, :, :].values, sal[90, :, :].values, iso)
Now, Let's plot it in Cartopy¶
In [47]:
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
lons, lats = np.meshgrid(lon, lat)
h1 = ax.contourf(lons, lats, isotemp, np.arange(0,16,0.5), trans=ccrs.PlateCarree(), cmap='RdBu',
vmin=2, vmax=14)
plt.colorbar(h1, orientation = 'vertical', fraction=0.02, label='ËšC').set_label(label='ËšC',
#weight='bold',
size=15)
ax.set_title('Temperature field along the \n 90 isopycnal mark', size=20);
Test Cell for Func 3: iso_average¶
In [37]:
temp_avr = iso_average(temp[90, :, :].values, temp_iso, ds.lev.values)
temp_avr
Out[37]:
In [48]:
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import AxesGrid
def make_map(projection=ccrs.PlateCarree(), figsize=(20, 15)):
fig, ax = plt.subplots(
figsize=figsize,
subplot_kw={"projection": projection})
return fig, ax
fig, ax = make_map(projection=ccrs.PlateCarree(), figsize=(20,15))
# fig, ax = plt.plot(projection = ccrs.PlateCarree(), figsize = (20,15))
# fig = plt.figure(figsize=(20,15))
# ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines(resolution="110m", color="k")
ax.add_feature(cfeature.LAND, facecolor="0.75")
lat = ds.lat.values
lon = ds.lon.values
#lons, lats = np.meshgrid(lon, lat)
h2 = ax.contourf(lon, lat, temp_avr, np.arange(0,20,2.0), extend='both', trans=ccrs.PlateCarree(),
cmap='RdBu')
plt.colorbar(h2, orientation = 'vertical', fraction=0.02).set_label(label='ËšC',
#weight='bold',
size=15)
ax.set_xticks(np.linspace(-180, 180, 5), ccrs.PlateCarree())
ax.set_yticks(np.linspace(-90, 90, 5), ccrs.PlateCarree())
ax.set_title('Vertically averaged temperature\nabove the 90 isopycnal mark', size=20);
# lon_formatter = LongitudeFormatter(zero_direction_label=True)
# lat_formatter = LatitudeFormatter()
# ax.xaxis.set_major_formatter(lon_formatter)
# ax.yaxis.set_major_formatter(lat_formatter)
In [ ]:
In [ ]: